function [AUC,VAL]=roc(obs,pred,plt,salto)

% ROC Response Operating Characteristic curves
% Syntax: AUC=roc(obs,pred,{plt},{salto})
% Description:
%    Creates a ROC curve plot of sensitivity vs 1-specificity
%    which is equivalent of p(true positives) vs p(false negatives)
%    
%    Input Arguments:
%
%    obs - observed outcome
%    pred - predicted outcome
%    plt:
%        1 - regular ROC curve plot
%        2 - ROC curve with threshold values
%        any other option will cause the graphic
%        display not to be produced.
%
%    Output arguments
%
%    AUC - area under the curve
%    VAL - values definign the ROC curve 
%
% Jonas Almeida, May 2006

if nargin<3;plt=0;end % By default do not plot
if nargin<4;salto=unique(pred)';end

obs_pos=sum(obs==1); % Observed positives = 100% true positives
positivos=pred(find(obs==1));npos=length(positivos);
negativos=pred(find(obs==0));nneg=length(negativos);

VAL=[]; 
for tr=salto    % tr= threshold
    VAL=[VAL;[tr,sum(negativos>=tr),sum(positivos>=tr)]];
end
VAL(:,2)=VAL(:,2)./nneg;
VAL(:,3)=VAL(:,3)./npos;
% inverting order
VAL=VAL(end:-1:1,:);

% Area under the curve:
arestas=[VAL(1:end-1,2),VAL(2:end,2),VAL(1:end-1,3),VAL(2:end,3)]; % edges
AUC=sum((arestas(:,2)-arestas(:,1)).*(arestas(:,3)+(arestas(:,4)-arestas(:,3))./2)); %area under the ROC curve

%Graphic display of ROC curve
switch plt
case 1
    h=area(VAL(:,2),VAL(:,3));
    set(h,'FaceColor',[0.75 0.75 0.75])
    hold on
    plot([0;1],[0;1],':');
    ylabel('sensitivity')
    xlabel('1-specificity')
    text(0.6,0.4,[' Area = ',num2str(AUC)],'FontSize',14)
    gridover
    hold off
    
case 2 % plot threshold value too
    subplot(2,1,1)
    h=area(VAL(:,2),VAL(:,3));hold on
    set(h,'FaceColor',[0.75 0.75 0.75])
    %set(h,'LineWidth',2)
    %plot(VAL(:,2),VAL(:,3),'-','LineWidth',2)
    ylabel('sensitivity (true positives)')
    xlabel('1-specificity (false positives)')
    hold on
    plot([0;1],[0;1],'--k','LineWidth',2,'Color',[0.2 0.2 0.2]);
    hold off
    axis=[0 1 0 1];
    subplot(2,1,2)
    plot(VAL(:,2),salto(end:-1:1)','LineWidth',3,'Color',[0.5 0.5 0.5])
    G=get(gcf);
    %set(G.Children(1),'XAxisLocation','top')
    ylabel('segmentation value (decision threshold)')
    G=get(gcf);G1=get(G.Children(1));G2=get(G.Children(2));
    set(G.Children(1),'Position',[0.1,0.1,0.8,0.8],'Color','none','YAxisLocation','right','XLim',[0 1]);set(G.Children(2),'Position',[0.1,0.1,0.8,0.8],'XLim',[0 1],'YLim',[0 1],'Color','w')
    YL=get(G.Children(1),'YLim');
    AUCstr=[num2str(AUC),'000'];
    text(0.5,YL(1)+0.05*(YL(2)-YL(1)),[' Area = ',AUCstr(1:4)],'FontSize',14,'Color','w')
    axes(G.Children(2));gridover;axes(G.Children(1))
end


% ----- call in functions -----
function gridover() % grid overlay
%axis([0 1 0 1]);
ax=[0 1 0 1];
G=get(gcf);G1=get(G.Children(1));hold on
for i=1:length(G1.XTick) % X grid
    plot([G1.XTick(i),G1.XTick(i)],[ax(3),ax(4)],'k:')
end
for i=1:length(G1.YTick) % X grid
    plot([ax(1),ax(2)],[G1.YTick(i),G1.YTick(i)],'k:')
end
hold off
        